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ABSTRACT 

This paper describes a cell averaging method for the ChebyBhev approximations of first 
order hyperbolic equations in conservation form. We present formulas for transforming 
between pointwise data at the collocation points and cell averaged quantities, and vice- 
versa. This step, trivial for the finite difference and Fourier methods, is nontrivial for the 
global polynomials used in spectral methods. We then prove that the cell averaging methods 
presented are stable for linear scalar hyperbolic equations and present numerical simulations 
of shock-density wave interaction using the new cell averaging Chebyshev methods. 
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1 Introduction 


In this paper we introduce a new spectral technique for the numerical solution of nonlinear 
hyperbolic equations. This technique, as almost every modern finite difference scheme for 
shock computations, is based on the cell averaged form of the equations. This is essential for 
finite difference shock capturing techniques and it is our experience that it plays an essential 
role in a successful spectral simulation of problems that involve shock waves [1]. 

Consider the nonlinear hyperbolic equation 

/ u, + F.(U) = 0, X € [-1,1] ,, n 

\ tf(x,0) = «,(.), ' 

with appropriate boundary conditions. 

The cell averaged form of (1.1) is obtained by integrating (1.1) between any two points 
— 1 < a < b < 1 to get 

&=mrr a J! ^ = -rraWM ~ (12) 

Let u(x,t ) be an approximation to at time t. Following Harten [4] we express the 

approximation to U(x,t + t) by 

u(x,t + r) = AE(t)H(-]u) y (1-3) 

where A is the cell averaging operator and E(t) is the exact time evolution operator corre- 
sponding to (1.1). Throughout the paper we will not distinguish between and 7 Zu. 

The operator 7Z(-; u) is of extreme importance, it represents the way we reconstruct u from its 
given cell average values u 7 _i = /*/ _l u{x)dx, where {xj}^ =l are the grid points. For 

finite difference schemes u is a piecewise polynomial of low degree, so that the reconstruction 

itself is simple. It becomes complicated only if one imposes also the requirement that the 
reconstruction should be essentially nonoscillatory. In [1] we have presented an essentially 
non-oscillatory Fourier method based on the cell averaging formulation (1.2). In that case 
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the transformations between the cell averages and the point values are simple and can be 
carried out efficiently by the Fast Fourier Transformations (FFT). This can be attributed to 
the fact that the boundary conditions are periodic and that the cell average of a trigonomet- 
ric function is proportional to the function itself . However, for Chebyshev methods, the cell 
averaging operation (denoted by the operator A) is not simple nor iB 7Z (- ; u ). As a matter of 
fact not only the formulation but also the implementation of A and 71 is not straightforward. 

In this paper we formulate the cell averaging technique for the Chebyshev method. We 
will discuss its stability for linear problems and show an example of its applicability to 
nonlinear systems of equations by simulating the problem of shock-density wave interaction. 
The cell averaging formulation is an essential part of the numerical code. 

The outline of the paper is as follows: In Section 2 we show how to reconstruct efficiently 
point values of a polynomial from its cell averages and vice versa.In Section 3 we introduce 
the new numerical technique and show its stability for linear problems. Section 4 is devoted 
to numerical results obtained by using the new method. 

2 Cell Averages and Point values 

In this section we will discuss the cell averaging operator A and the reconstruction operator 
71 in the context of the Chebyshev methods. In these methods the approximations are 
taken from the space of polynomials of degree N. It is therefore clear that A and 71, when 
restricted to the polynomial space, can be expressed as matrices An and Rn. We will give 
the explicit formulas for these matrices. We start by discussing the operator A. 

Assume that /(x) is in C r [— 1, 1], r > 0. Let Xj = cos(^), 0 < j < N be the Chebyshev- 
Lobatto points in [-1,1]. For later use, define x 7 _i = cos(^j^), 1 < j < N. The cell 
averaged function /(x) of f(x) is defined as 

1 Ma(ss) 

J(x) = Af =: -r—r-T r-7 — r- / f(x)dx for -1 < x < 1, (2.1) 

h 2 (x) - hi{x) Jh^x) v J 
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where 


/ii(x) = cos (cos x — J, 

, , , , A0, 

h 2 [x) = cos (cos x + — -), 

it 

. _ 7T 

A0 = — . 

iV 


( 2 . 2 ) 


The reason for the definition in (2.1) is that 

/(*,•_ i) = / ; /(®)dx, for 1 <j<N. (2.3) 

a Xj_! — Xj Jxj 

As stated before, we are interested in A operating in a polynomial space. In Lemma 1, 
we show that the result of A on a polynomial is still a polynomial of the same degree. 


Lemma 1 Let U k {x ) = i( x )> k ^ 0, be the second kind of the Chebyshev polynomials. 


Then 


where 


V k {x) = AU k = cr k U k (x), 


°k = 


sin(& + l)f 

(* + i) sin ¥ 


(2.4) 


(2.5) 


Proof: Substituting U k (x) in the right hand side of (2.1) and making the transformation 
x = cos 6, 0 < 6 < 7r, we have 

1 yco«(co« -1 x+^p) 

^ k ^ ^ cos (cos -1 x + f) — cos (cos -1 x — f) J c<M ( C ot - 1 x-^-) ^ 

Since T k (x ) = cos(kd), then U k (x ) = and therefore 


dx. 


U k (x) = 


r e +7 

/ ^ sin (A: + 1)5 d6 

JS—Qr 


2 (k + l)sinf sin0 

1 


cos ( k + 1)0 


2(k + l)sin f sin0^ ^ A: + 1 


e+¥ 


sin ( k + 


( k + l)sin 


l)f ( 1 sin (A 

inf \k+ 7 ri 


(A + 1)0 N 

sin 6 
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Q.E.D. 


i.e. 

U k (x) = a k U k {x). 


^From Lemma 1 one gets 

Corollary 1 The cell averaged function of any polynomial of order N is a polynomial of the 
same order. 


Proof: any polynomial of degree N has the expression 

/(*) = £<■»£«*), (2.6) 

k=0 

where a* are constants. 

Therefore, by Lemma 1 

f(x) = X) ° k a-kUk{x). (2.7) 

k = 0 

Q.E.D. 

Thus, Lemma 1 gives explicitly the eigenvalues cr* of the matrix (the restriction of 
A to polynomials of degree N ). These eigenvalues are uniformly bounded from above and 
below, in fact 

- < a k < 0 < k < N. (2.8) 

7T 2 

If f(x) is a polynomial of degree N-, then it is uniquely determined by its values f(xj), j = 
0, • • • , N. So theoretically f{xj_i), j = 1, • • • , N can be determined. Therefore the transfor- 
mation from f(xj), j = 0, • • ■ , N to /(a^ i ), j = 1, • ■ ■ , N is well defined and we only need 
to address the issue of its efficient implementation. 

In general it is known that 

/(*)=E^(x), (2-9) 

k=0 
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where 


^ = ■=£■£; /mum. 

j=o 

Co = Ci\r = 2, c*. = 1 if A: ^ 0 or N. (2.10) 

Alternatively 

/(*) = £ W*), (2.H) 

Jfc=0 

b k = k = N -l or W, 

z 

6* = ^(cfcO* - ajt+a), 0 < k < N -2 (2.12) 

and therefore by corollary 1 

/(*) = X>* Ws). (2- 13 ) 

jfc=0 

Now {/(®j_^)}jLi are obtained by substituting m (2.13) (this can be carried out 
using the FFT). 

We note that equations (2.9) - (2.13) describe how to get the vector /(»,■_ i), j = 1, • • ■ , AT 
from the vector /(x<), t = 0, • • • , N. Denote by A N the Nx(N + l) matrix defined by this 
transformation. We note that An can be written explicitly. In fact the polynomial /(x) has 
an unique representation as 

/(*) = £ f( x j)9j{ x )> ( 2 - 14 ) 

3 = 0 

where 

, , (1 - * ! )TJ,(*)(-1)> +1 2 ", r t (»j)T t (») 

- ij) 5 ,-at 

with 2^ defined in (2.10) . 

It follows upon substituting (2.14) in (2.1) and using the fact that f 0 (x) = U 0 {x), f\(x) = 
2 a x Ux{x) and T k (x) = \{o k U k (x) - <J k ^U k . 2 (x)) for k > 2, 

I{ x ) = T,K x t)9l{ x )> ( 2 - 15 ) 

i=o 
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where 


9i{ x ) — 2 + 2<TiTi(x/)Ui(x) + 


2 ^ n(^)(^(x)-a fc _ 2 ^_ 2 (x)) 

fyN k=2 2c k 


with cr* defined in (2.5). 

Setting x = Xj_i in (2.16) 


(2.16) 


(An)* = i), 0 < i < N, 1 < j < N. (2.17) 

Thus we have outlined two procedures to get J{xj_ i) from /(xj), one uses the FFT and 
another uses matrix vector multiplications. 

We are now ready to discuss the reconstruction operator 7Z (- ; /). Note that in the 
beginning of the solution process (1.3), we only have the values fj_i , j = !,-■ ■ } N, thus we 
need another piece of information in order to define uniquely the N - th degree polynomial 
/(*)• This piece of information is provided by the boundary condition. For simplicity, we 
assume that the boundary condition is of the form 


/(!) = /(x 0 ) = / 0 . 


(2.18) 


The reconstruction is done in two steps. Define first a (N - l)-th polynomial /jv-i(s) 
which collocates f(x) at , i.e. ftf-i{xj_i) = /(x^_ i), 1 < j < N, it is readily 

verified that 


JV-l 

/jV-l(®) = c kTk(x) , 

k=0 


where 


N 


c k = -z^Ylf( x j-L)Tk(xj_ i), c k is same as in (2.10). 
Alternatively we have 


/tf-i(x) = £■ b k U k (x), 

k=0 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


where 


h =j, k = N- 2, JV-l, 

h = ^(cfcc fc - c fc+2 ), 0<A;<JV — 3. (2.22) 
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Now, by Corollary 1 


(2.23) 


/jv-i(x) = 53 

k= o a * 

is a polynomial of degree N — 1 such that Af^-i = fn- i- 

Generally, /n-i( x ) does not satisfy the boundary condition (2.18). There are two ways 
to modify /jv_i(s) so that the boundary condition (2.18) is satisfied. In the first way we can 
add to /n_i(x) an N - th degree polynomial Q(x) such that 


Q (* j -0 = °, !<;<*, 


and 


f N-l ( x ) + Q(l) — fo- 

It can be verified that in order to satisfy (2.24) 

<?(*) = c(( l-x J )Ti(x))'. 

Let f(x ) be the sum of /at_i(x) and Q(x), 

/(*) = /*.,(*) + c ((1 - *’)!*„(*))' 

= - c[xTi,(x) + N 1 T„(x)]. 


(2.24) 


(2.25) 


(2.26) 


(2.27) 


The last equality follows from the Chebyshev equation, and the constant c is now deter- 
mined by the condition /(x o) = fo, i-e. 

1 


c — — 


2 N 2 


(fo ~ /n- i ( 1 ))- 


(2.28) 


Finally given f(xj_ i ), j = 1 , • • • , N and fo we can get 

/(*i) = /*-,(*;) - c , i = 0,--,N (2.29) 

where /jv_i(xj) can be evaluated from (2.23) by using the FFT. 

Note that in this procedure we change the values of fN-i( x i) a t all the grid points. 


7 


Denote by Rn the ( N + 1) x ( N + 1) matrix transforming f 0 and f(xj_ i), 1 < j < N 
to 0 < j < N, i.e. 

ifM,-” ,f(x N )) J = R n (/(*o),/(*j),*-*,/(*jy_j)) T - (2.30) 

As before we can write Rn explicitly. Equation (2.19) can be rewritten as 

In- i(») = (2.31) 

i=i 

where hj(x k _±) = Sjk and Aj*(z) are polynomials of degree N — 1; explicitly 

s j( x) = (-iy»i V .in(y-I)^)^a r . (2.32) 

3 ~ 2 

It can be shown that ^j(x) is the cell averaged polynomial of 

hj(x) = £ A h U k (x) (2.33) 

k=o 

with A* defined as 


A* = if fc = JV-2.JV- 1 , 

At = j^(3t(x ) _i)-It + j(x i _i)), ifO<t<JV-3. 

As a result of (2.31) and (2.33) polynomial /jv_i(x) takes the form 


In- i(x) = J2 f(xj_i)hj(x), 

J-l 


(2.35) 


/(x) = /w-i(x) - ± (/„ - /»-,(!)) ((1 - *’)?”»(*))' 


Using (2.31) we get 


/(*) = - Tm ( /o - £/(*<-»)*»•(!)) ft 1 - (2.36) 

= Ef( x i-i) *iW + ^(( 1 - I ’) I »( 1 ))' - 2^ ((1 - x’Wt*))' 
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Substituting x = x< into (2.36) gives 

«i i, if i = 0, 

-ji, ((1 - x? )T K ( Xi ))' , if 3 = 0, 

M * i ) + |#((1 - *?)!«*<))'. 1 < f,J < 

To summarize we state 


(Rn) 0 - = 


(2.37) 


Lemma 2 Let /(x) = ft(- ;/) 6e tfie JV-t/i polynomial defined in (2.27) then 

(^4/)(x j _i) = f{x j _ i) fori <j <N. (2.38) 


A different way to modify /jv_ i(x) in (2.23), in order to satisfy the boundary condition 
(2.18), is to add to it a polynomial Q x (x) of degree N such that the point values /y_ x (x,) 
remain unchanged and the new polynomial satisfies the boundary condition. Thus instead 
of (2.27) we define 

f(x) = f N - i(x) + (/o - /jv-i(I)/ + 2J\P^ ( 2 ‘ 39 ) 

Computationally (2.39) is simpler than (2.23) (2.27). However, note that in this case 

( 2 - 4 °) 

which is in contrast to (2.38). The matrix corresponding to (2.40) can be formed similarly 
as in (2.37). 


3 Cell Averaging Chebyshev (CAC) Method and 
Linear Stability 

In this section we will establish the stability of the Cell Averaging Methods, presented in 
Section 2, when applied to a first order scalar hyperbolic equation. It is tempting to try to 
obtain stability in the L 1 norm because of the way the method (1.3) is presented. However 
we will only give the stability estimate based on a weighted Chebyshev norm. 
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Consider the initial boundary value problem of the scalar hyperbolic equation 

r U t = U x , *€[-1,1], 

< U(x, 0 ) = U 0 (x), ( 3 . 1 ) 

l u(i,t) = o. 

The cell averaged form of (3.1) is 

0 + h,( x) - u{K{ X ), o] = o, (3.2) 

where h x (x) and /^(x) are defined in (2.2). 

We follow the notation (1.3). Let the jV-th degree polynomial u(x,t ) be the approxima- 
tion to (3.2). Our aim is to find the error equation for TZ(-;u) defined either by (2.27) or 
(2.39). iFrom (1.3) 

u(x, t + r) = AETZu(x, t). (3.3) 

Applying the reconstruction operator TZ we get 


7Zu(x,t + r) = TZAETZu(x,t) 


(3.4) 


iProm equation (3.4) it is clear that TZu satisfies exactly the equation 


dTZu 

~df 


dTZu 

dx 


+ r [(!-*%,]' 


for the reconstruction procedure (2.27) and 


dTZu dTZu 
dt dx 


+ 7l(l + x )T'n 


(3.5) 


(3.6) 


for (2.39). t and t x are quantities depending on time t. 

Note that (3.6) is the error equation for the Chebyshev collocation method. For a scalar 
linear equation the C AC method corresponding to (3.6) is equivalent to the known collocation 
method [2], It remains to investigate the stability of (3.5). For simplicity, in the remainder 
of this section we denote T^u by u. From the construction of 7Z in section 2 we know that 
u(x } t) satisfies the boundary condition in (3.1), i.e. u(l,<) = 0. 
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It is interesting to note from (3.6) that the CAC method with reconstruction operator 
(2.27) is equivalent - for constant coefficient case - to the collocation method where the grid 
points are the zeros of the polynomial 


<?(*)= [(i-x J )r w (*)]'. 


Note that by using standard identities 


Q(x) = - [xT^(x) + JV 2 7V(x)] . (3.8) 

The term r in the right hand side of (3.5) is determined by the boundary condition 
u(l,t) = 7£ti(l,t). In fact substituting' x = 1 in the equation (3.5) and noting that |^(1, t) = 


0 we get for r 


«.(M) 


An alternative expression can be obtained by equating the highest coefficient on both 
sides of (3.5), thus if u is expanded as 


u(x) = u k T k (x) 


(3.10) 


_ 1 dupr / 

T “ ~N(N + 1) dt 1 

Before stating the main stability result of this paper, we state the following lemma 


(3.11) 


Lemma 3 If f(x) = Y^t=o 1 a kT k {x), then 




where Cj is defined in (2.10); Xj 0 < j < N are the Chebyshev-Lobatto points. 


Proof: see [3]. 

The main stability result follows from the following lemma 
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Lemma 4 Let U(x,t) be the solution of (3.2) and u{x,t ) = IZu be the approximation to 
U(x,t) by the CAC method (1.3) with the reconstruction operator ( 2.27) or (2.30), then 


d_ 

dt 


“X 1 Xj 2 f j\ i 7r -^ «2 f±\ 

— L ;Ti ~ U ( x i» 0 + o ? xr , ix MO 


2iV Cj 1 - Xj 


2(N -f 1) 


< 0. 


(3-13) 


Proof: We multiply (3.5) by i) and sum from the points xn = —1 to x x to get 


8 7T ^ 1 1 + Xj 2 


7T ^ 1 1 + Xj 


V n y— v X J. T «L 7 9/ v I 1 T *1 / \ / 


- (3.14) 

iV i=! c i 1 ~ x i 


We treat the two terms in the right hand side of (3.14) separately. 


T 7T y ^ 1 1 “|“ (Z)« . . , , 

7 = — —“(*>> *K(*j»*) 

iV J=1 c i 1 “ 

7T y — -s 1 1 ~j” (C* , , . 7T b. . 

= jy 1 _ a ^ u ( a b> 0 u *( z j> *) + jy u x(M)- 
By noting the exactness of the Gauss Lobatto formula one gets 

y f 1 1 *f X . . 6 . . dx 7T . 

1 = L — x u{x ’ () 7r^P + jv“- (1 - 1 

and integration by part yields 

r 1 i + f ^ 2 (x ,0 


r _ _ j + 1 , , *, n ,x 

7 - vT^(l - x) 2 ^ + N^ 1 ’^ 


We use again the Gauss Lobatto formula to get 


7 *£*< 1+ 2V xj 2 + JV^ 1, ^ 


7/ 


X 


JX tt 2 (gj»*) , » a/ 


and therefore 




We turn now to the second term in (3.14) 


II = -rJV’£ £ |-±±fi u{ly , <)T„(X,) 

iV j=i c j 1 *j 


(3.15) 


(3.16) 


(3.17) 


(3.18) 


(3.19) 


(3.20) 


12 



or 


II = -Tl i - rW J ^,(l,t). 
™ j=o c i 1 


(3.21) 


Using (3.9) and (3.11) for r one gets 

N du N 7T ^ 1 1 + X j _ JL„V 1 *Y (3.22) 




7T 1 1 H” Xj ✓ i\m *.2/1 

T-af-ArE^i -fuix^t^Xj) - 2^(1, t). 

1 at N Cj 1 - Xj ZiV 

As I±|u(x,t)Ttf(s) is a polynomial of degree less than 4 N — 1, one gets using lemma 3 
ji = £ i + 7ra,N ' (3 ' 23) 


where anf is the 2 N coefficient in the expansion of t)7Vr(®) 

It can be easily verified that 




02AT = — 


BO 


7/ = — 7T 

Using (3.19) and (3.26) we get 


JNT . du* it 2 , v 

FTT UAr "9r " 2iv x( ’ } ‘ 


iV dtlff 7T j/. ,\ 
I + II < -7TT7— T^N— - Tat U *( 1 > 0- 


(3.24) 

(3.25) 

(3.26) 

(3.27) 

N + l~" dt AN~ X ''~’ ~ r 
Substituting (3.27) into (3.14) yields (3.13). This completes the proof of the lemma. 

Q.E.D. 

We note that from Lemma 3 

h t () " 5 /-. r lm - (328) 

so that we can finally state that 

Theorem 1 Let u(x,t ) = Hu be the solution of the GAC method with the reconstruction 
operator (2.27) or (2. SO), then 


1 r 1 1 + x u 2 (x,t ) 


1 r l± 

2 7-i 1 — 


s/l — x‘ 


dx 


27V — 1 
4(JV + 1) 


2 m< I f 1 1 + a? u 2 (a,0) dx 

~ 2 J- 1 1 — x s/1 — x 2 4 (N + 1) 


2N — 1 


7TU^(0). 

(3.29) 
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4 Numerical Results 


In this section we apply the CAC spectral method (1.3) to the one dimensional gas dynamics 
equations. The time evolution is done by the Runge-Kutta type method. Each step of the 
Runge-Kutta scheme is done as follows: 


Fully Discretized CAC method 
Step 1 : Reconstruction: 

given u"_i , j = 1, • • ■ , N } we use the boundary condition and the matrix Rn to find the 
point values uj, j = 1, • • • N as suggested in (2.27) or (2.30); 

Step 2 : Solution in time: 

we update the values j = 1, • • ■ , N using the forward scheme, 

J y 

At 


= a U - zrb- ( F ^> - f <“?)] ,i = 


(4.1) 


where At is the time step. 

The reconstruction operator Rjy yields spectrally accurate point value approximations 
to the exact solutions if the exact solutions are smooth, thus Uj_ i, j = can be 

expected to approximate their cell averages accurately. However, if the exact solutions are 
discontinuous, the point values Uj , j = 0, • • • , N obtained by Rn will be oscillatory as the 
result of the Gibbs phenomenon. In [1] we proposed a practical way to obtain essentially 
nonoscillatory spectral reconstruction to a discontinuous function from its oscillatory Fourier 
approximations. The key idea there is to augment the Fourier space by adding simple dis- 
continuous functions whose locations and magnitudes of discontinuities are approximations 
to those of the shock waves in the numerical solutions. In our computations of CAC method, 
we extend this idea along the same line to obtain essentially nonoscillatory reconstructions. 
The estimates on these reconstructions will be appearing in a separate work. We refer the 
reader to ([1]) for further details. 

Now consider the Riemann problem for the Euler equations for a polytropic gas 
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Ut(x,t) + f*(U) = 0, 


(4.2) 


where U(x, t) and f(U) are defined as 

U = (p,M,E) T , -1 <x < 1, 

f(U) = qV + (0, P,qP) T , (4.3) 

where 

P = (7 - 1)(£ - JM 1 ). M = pq, (4.4) 

with 7 = 1.4 and the initial conditions are as follows 

( pL,qL,PL ) = (3.857143,2.629369,10.33333) when x < -0.8, 

( PR,qR,PR ) = (1 + esin57rs,0,l) when x > -0.8, (4.5) 

where e = 0.2. 

The solutions to (4.2) - (4.5) model the interaction between a moving shock wave and 
disturbances. Note that in the right state of the density a sinuous perturbation of magnitude 
e = 0.2 is superimposed upon a constant state. ^From linear analysis it can be shown that 
the disturbances will interact with the shock wave. A density wave of different frequency 
and magnitude will emerge behind the shock wave. Also the disturbance in the density 
field will perturb the velocity and pressure fields behind the shock wave. The numerical 
solution of this Riemann problem mandates a high order scheme in order to capture the fine 
structures in the solutions for a correct interpretation of the physical process. We test this 
problem with second order MUSCL scheme [5] and third order point value version ENO finite 
difference scheme [6] and the CAC spectral method proposed in this paper. Our numerical 
results have shown clearly the advantage of a higher order numerical scheme for problems 
with complicated structures. 

Figure 1-3 show the density profiles obtained by the three methods mentioned above. 
Figure 1 is the result using CAC spectral method. The solid lines are the solutions taken 
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from [6] using the third order ENO finite difference method with N = 800 which we take as 
the exact solutions . Figure 2 and Figure 3 are the results with the second order MUSCL 
scheme [5] and the third order ENO finite difference scheme respectively. In all three cases 
we use the same amount of mesh points N = 220. All the results are plotted at the same 
time t = 0.3. 
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Figure 1: The CAC spectral method: density, N = 220, time t = 0.3. (+) - numerical 
solutions, solid line - exact solutions. 
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Figure 3: The third order ENO scheme: density, N = 220, time t = 0.3. (+) - numerical 
solutions, solid line - exact solutions. 
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